R/Function 9 oceanVulnerability v.0.04.r

In monteroserra/rOceansNew: Spatial Trends in Marine Biodiversity and Conservation

#' Function 9 oceanVulnerab: link hotspots of biodiversity to projected climatic impacts
#'
#' @description access and links climatic projections to hotspots of biodiversity to unravel potential impacts of climate change 
#' @param biodiversity_grid raster layer with any biodiversity metric previously computed with oceanDiversity
#' @param biodiversity_metric the name of the biodiversity metric
#' @param reshape_climatic_layer whether climatic layer resolution must be changed to match biodiversity grid resolution
#' @param new_cell_size define desired cell size of the climatic layer for the analysis
#' @param map_climtatic_impacts whether we want the resulting spatial grid of future trends to be mapped
#' @param plot_histograms whether histograms of expectec temperature increases at hotspots are to be mapped
#' @return Object of hotspots coordinates and climatic projections according to IPCC scenarios RCP 2.5, RCP 4.5, and RCP 8.5
#' @details This function allows accessing and exploring expected changes environmental parameters 
#' @export


oceanVulnerab = function (biodiversity_grid,
                          climatic_layer =NULL,
                          get_climatic_layers = T,
                          reshape_climatic_layer = F,
                          new_cell_size = 1,
                          plot_histograms = T,
                          map_climtatic_impacts = T,
                          low_color="steelblue", 
                          mid_color="gold",
                          high_color= "firebrick", 
                          col_steps=20) {
  
  if(reshape_climatic_layer){
    
    factor_conver = new_cell_size / res(climatic_layer)[1]
    climatic_layer <- aggregate(climatic_layer, fact=factor_conver)
    
  }
  
  
  hotspots = as(as(biodiversity_grid, "SpatialPointsDataFrame"), "data.frame") 
  
  colnames(hotspots) = c("hotspot", "x", "y")
  
  xy = SpatialPoints(hotspots[,c("x","y")])
  projection(xy) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
  
  if(get_climatic_layers){
    
    climatic_layerRCP8.5 = oceanFuture(IPCC_scenario = "RCP85")
    climatic_layerRCP4.5 = oceanFuture(IPCC_scenario = "RCP45")
    climatic_layerRCP2.5 = oceanFuture(IPCC_scenario = "RCP26")
    
    hotspots$future_SST_RCP8.5 = extract(climatic_layerRCP8.5, xy)
    hotspots$future_SST_RCP4.5 = extract(climatic_layerRCP4.5, xy)
    hotspots$future_SST_RCP2.5 = extract(climatic_layerRCP2.5, xy)
    
    
    if(plot_histograms){
      
      par(mfrow=c(1,3)) 
      
      h25 = hist(hotspots$future_SST_RCP2.5, breaks=15, main="RCP 2.5",col="steelblue",
                 xlab="Expected Increases in Tª(ºC)", ylab="Hotspots of Biodiversity")
      box()
      
      h45 = hist(hotspots$future_SST_RCP4.5, breaks=15, main="RCP 4.5",col="orange",
                 xlab="Expected Increases in Tª(ºC)", ylab="Hotspots of Biodiversity")
      box()
      
      h85 = hist(hotspots$future_SST_RCP8.5, breaks=15, main="RCP 8.5",col="firebrick",
                 xlab="Expected Increases in Tª(ºC)", ylab="Hotspots of Biodiversity")
      box()
      par(mfrow=c(1,1))
    }
    
  } else {
    
    hotspots$future_layer = extract(climatic_layer, xy)
    
    hist(hotspots$future_layer, breaks=15, col="firebrick",
         xlab="Expected Changes", ylab="Hotspots of Biodiversity")
    
  }
  
  
  
  if(map_climtatic_impacts){
    
    dev.off()
    
    color.gradient <- function(x, colors=c(low_color,mid_color,high_color), colsteps=col_steps) {
      return(colorRampPalette(colors) (colsteps) [ findInterval(x, seq(min(x),max(x), length.out=colsteps)) ] )
    }
    
    palette(color.gradient(1:20))
    
    plot(biodiversity_grid, col="white", legend=F, main="Climatic impacts RCP8.5")
    maps::map("world",add=T, fill = T,bg="grey20",col="grey30")
    points(hotspots[,c("x","y")], pch=21,cex=1.5, 
           bg=factor(hotspots$future_SST_RCP8.5))
    
    
  }
  
  return(hotspots)
  
}